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FOREWORD 

Thi^  report  examines  multiple  satellite  orbit  determination  for 
neighboring  satellites.  The  procedures  discussed  can  result  in  a 
reduction  of  both  the  total  number  of  parameters  and  the  computation 
time  required  to  process  satellite  ephemerides  individually. 

The  work  described  in  this  report  was  performed  in  the  Astronautics 
and  Geodesy  Division,  Warfare  Analysis  Department.  The  research  was 
funded  by  the  Naval  Electronics  System  Command  under  SPACETASK  PME-106- 
511-009C5X3519. 

The  report  was  reviewed  by  R.  J.  Anderle,  Head,  Astronautics  and 
Geodesy  Division. 


RELEASED  BY: 

y y/.  a ) u ‘ 

R.  A.  NIEMANN,  Head 
Warfare  Analysis  Department 


UINLJjttSOir  i£U 


/ 

1 


SECURITY  CLASSIFICATION  OF  THIS  PAGE  (Whin  Data  Enlarad; 


i-  REPORT  DOCUMENTATION  PAGE 

READ  INSTRUCTIONS 
BEFORE  COMPLETING  FORM 

1 REPORT  NUMBER 

TR-3308 

2.  GOVT  ACCESSION  NO. 

1.  RECIPIENT'S  CATALOG  NUMBER 

«.  TITLE  (»nd  Submit) 

MULTIPLE  SATELLITE  ORBIT  DETERMINATION  FOR 
NEIGHBORING  SATELLITES 

s.  type  of  report  a period  covered 

S.  PERFORMING  org.  report  NUMBER 

7 author^*; 
Patrick  J.  Fell 

i.  contract  or  grant  number^; 

9.  PERFORMING  ORGANIZATION  NAME  ANO  ADDRESS 

Naval  Surface  Weapons  Center 
Dahlgren  Laboratory 
Dahlgren,  Va.  22448 

to.  PROGRAM  ELEMENT.  PROJECT.  TASK 
AREA  & WORK  UNIT  NUMBERS 

1 1 CONTROLLING  OFFICE  NAME  AND  ADDRESS 

12.  REPO  r D a * E 

April  1975 

13.  NUMBER  OF  PAGES 

1 14  MONITORING  AGENCY  NAME  ft  ADDRESS  (It  different  from  Controlling  Office) 

15.  SECURITY  CLASS,  (ol  Ihlt  "port) 

UNCLASSIFIED 

15a  DECLASSl  F|  CATION;  DOWN  GRADING 
SCHEDULE 

16.  distribution  statement  roi  this  r»p art) 


improved  for  Public  Release;  Distribution  Unlimited 


17.  DISTRIBUTION  STATEMENT  (of  the  abatrect  entered  In  Block  20,  It  different  from  Report) 


18.  SUPPLEMENTARY  NOTES 


19  KEY  WORDS  (Continue  on  raver te  elde  If  neceeeery  end  Identify  by  block  number) 


20.  ABSTRACT  (Continue  on  reveree  elde  It  neceeeery  end  Identify  by  block  number) 

If  satellites  follow  orbit  paths  within  one  hundred  kilometers  of  each 
other, 

(i)  The  simultaneous  processing  of  orbital  data  using  a multiple 
satellite  orbit  determination  method  produces  ephemerides  for  the  satellites 
with  an  accuracy  consistent  with  individually  determined  orbits.  Hcwever 
this  method  predicts  the  relative  separation  of  the  satellites  more  accurately 
in  the  presence  of  modeling  errors.  PRICES  c'  l ECT  TO 

DD  i JANM?3  1473  edition  OF  1 NOV  65  IS  OBSOLETE  UNCLASSIFIED 

S/N  0 1 0 2*  0 14"  660  1 • 


I 


SECURITY  CLASSIFICATION  OF  THIS  PAGE  (Whin  Data  tntind) 


Ul^Cl^DDXr  l£iU 


V.IJWITV  CLASSIFICATION  OF  THIS  *AGC(Wh#n  DM  Knffd) 


(ii)  The  use  of  perturbation  solutions  from  the  integration  of  one 
satellite's  variational  equations  to  process  the  other  satellite's  tracking 
data  introduces  a small  additional  error  into  orbit  solutions.  But  such 
a procedure  can  substantially  reduce  oarputation  time  required  to  ccnpute 
the  orbits  of  the  satellites. 

Based  on  this  analysis  it  is  reootmended  that  a multiple  satellite  data 
processing  scheme  be  adopted  for  satellites  in  neighboring  orbits.  In 
the  interest  of  accuracy,  hewever,  the  use  of  perturbation  solution 
approximations  should  be  reserved  unless  ccrrputation  time  requirements 
become  critical . 


SECURITY  CLASSIFICATION  OF  THIS  FAOEfWl*"  Data  Knit  tad) 


ABSTRACT 


If  satellites  follow  orbit  paths  within  one  hundred  kilometers  of 

each  other, 

the  simultaneous  processing  of  orbital  data  using  a multiple 
satellite  orbit  determination  method  produces  ephemerides  for  the 
satellites  with  an  accuracy  consistent  with  individually  determined 
orbits.  However  this  method  predicts  the  relative  separation  of  the 
satellites  more  aecurato.lv  in  the  presence  of  modeling  errors. 

fii)  The  use  of  perturbation  solutions  from  the  integration  of  one 
satellite's  variational  equations  to  process  the  other  satellite's 
tracking  data  introduces  a small  additional  error  into  orbit  solutions. 
But  such  a procedure  can  substantially  reduce  computation  time  required 
to  compute  the  orbits  of  the  satellites. 

Based  on  this  analysis  it  is  recommended  that  a multiple  satellite 
data  processing  scheme  be  adopted  for  satellites  in  neighboring  orbits. 
In  the  interest  of  accuracy,  however,  the  use  of  perturbation  solution 
approximations  should  be  reserved  unless  computation  time  requirements 
become  critical. 
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I.  INTRODUCTION 


Orbit  determination  for  a system  of  multiple  satellites  involves  the 
determination  of  an  extended  number  of  parameters,  usually  by  least 
squares  differential  correction.  In  some  instances,  depending  on  the 
orbital  geometry  of  the  system,  this  parameter  set  may  be  redundant.  In 
these  situations  parameters  that  are  common  among  various  satellites  may 
be  combined  to  reduce  the  total  number  of  state  parameters 
required  if  all  satellites  were  to  be  processed  individually.  These 
parameters  may  include  dynamical  quantities  as  well  as  bias  parameters. 

A least  squares  method  describing  this  parameter  reduction  approach  is 
given  for  satellites  in  neighboring  orbits.  The  method  is  then 
numerically  applied  to  satellites  which  are  100  kilometers  apart. 

Another  aspect  of  the  orbit  determination  problem  is  the  integration 
of  variational  equations  to  obtain  perturbation  solutions  necessary  for 
orbit  improvement.  These  solutions  are  normally  produced  by  a time 
consuming  numerical  integration  which  must  be  exercised  for  each  satel- 
lite in  the  system.  If,  however,  the  path  of  two  satellites  are 
close  to  each  other  then  it  may  be  possible  to  utilize  the  perturbation 
solutions  from  one  integration  for  both  satellites.  This  would  result 
in  a reduction  of  the  computation  time  necessary  to  process  multiple 
satellites.  An  examination  of  this  concept  is  presented  below. 


II.  Computational  Methods 


A.  Parameter  Reduction  Approach 

Consider  a system  of  multiple  satellites  orbiting  in  close 
proximity  with  a physical  separation  remaining  less  than  one  hundred 
kilometers.  With  this  particular  constellation  geometry  the 
satellites  are  viewed  simultaneously  by  each  tracking  station.  Under 
this  condition  a multiple  satellite  orbit  determination  scheme  can  be 
devised  based  on  the  assumptions  that  (1)  the  scaling  constant  for  each 
satellite's  atmospheric  <rag  model  are  equal  and  (2)  the  refractive 
structure  of  the  troposphere  is  equivalent  for  all  satellites  per  pass. 
These  assumptions  enable  the  total  number  of  state  parameters  to  be  de- 
creased when  all  satellite  data  are  processed  simultaneously.  This  re- 
sults in  a single  drag  scaling  parameter  for  the  constellation  and  only 
one  tropospheric  refraction  scaling  parameter  per  pass.  This  reduction 
in  the  total  number  of  estimated  parameter  results  in  a stronger 
covariance  for  the  computed  satellite  positions. 

The  formulation  presented  in  Appendix  B is  that  designed  for 
satellites  in  neighboring  orbits.  The  technique  presented  may  be 
generalized  in  an  obvious  way  for  other  satellite  systems  where 
parameter  reduction  is  possible. 


B . Perturbation  Solution  Approximation 


In  order  to  form  the  least  squares  normal  equations  en  route  to 
determining  the  best  ephemerls  for  a satellite  it  is  necessary  to  pro- 
duce certain  quantities  called  perturbation  solutions.  These  quantities 
as  mentioned  previously  arc  produced  by  a time  consuming  numerical 
integration  of  the  satellite's  variational  equations.  In  short,  the 
variational  equations  express  how  accelerations  experienced  by  the 
satellite's  center  of  mass  would  change  given  changes  in  orbital  con- 
stants at  epoch.  Numerical  integration  then  produces  a history  of 
variations  in  position  and  velocity  resulting  from  changes  in  initial 
orbital  conditions.  These  partial  derivatives  are  required  for  the 
application  of  the  chain  rule  in  computing  the  partial  derivative  of 
data  taken  at  time  t with  respect  to  a given  state  parameter  at  epoch. 

If  the  equations  of  motion  for  two  satellites  approximate  each 
other  then  the  differences  in  the  perturbation  solutions  for  the  two 
will  be  small  depending  on  how  "good"  the  approximation  is.  Therefore 
under  some  conditions  it  is  possible  to  utilize  the  perturbation 
solutions  from  one  integration  for  both  satellites: 

Let  tQ  be  some  epoch  time  and  let  { e j ( tQ) } ^ and  { ej ( tQ) ] ^ represent 
the  orbital  elements  for  two  satellites.  Assuming  that  the  satellites 
are  in  approximately  the  same  orbit  the  idea  is  to  determine  two  times 
t^  and  t^  when  the  quantity 

IVV  " T2(t2)l 


reaches  a minimum  during  the  first  revolution  after  epoch.  The  vectors 
r^  and  ^ represent  the  position  vectors  for  the  two  satellites.  Then 
using  numerical  integration  update  {ej(tQ)}^  from  tQ  to  t^  and  { ( tQ) } 2 
from  t to  t . The  time  t,  and  t„  become  the  new  epochs  for  the 

02  12 

respective  satellites. 

Now,  using  a given  step  size  At,  a perturbed  trajectory,  position 
and  partial  derivatives  versus  time,  is  integrated  for  satellite  1 for 
a desired  length  of  time.  An  inertial  trajectory,  position  and  velocity 
versus  time,  is  integrated  for  satellite  2 over  the  same  time  span.  A 
"synthetic"  perturbed  trajectory  is  then  created  for  satellite  2 by 
merging  its  inertial  trajectory  with  the  perturbed  trajectory  for 
satellite  1.  This  is  accomplished  by  extracting  the  time  (t£  + nAt)  and 
position  for  satellite  2 from  the  inertial  trajectory  and  the  partial 
derivatives  at  time  ft ^ + nAt)  from  the  perturbed  trajectory.  This 
"synthetic"  perturbed  trajectory  for  satellite  2 may  now  be  used  in  the 
least  squares  orbit  determination  process. 

Obviously  certain  error  levels  will  be  introduced  into  orbit 
solutions  due  to  this  approximation,  but  the  magnitude  of  the  errors 
depends  on  the  forces  acting  on  each  satellite  from  times  t^  and  t2  on. 
Numerical  examples  are  discussed  below  and  error  estimates  for  this 
technique  are  given. 
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III.  Numerical  Examples 


The  computational  procedures  described  above  were  applied  to  a 
system  of  three  coplanar  satellites  in  quite  similiar  orbits.  The 
satellites  maintain  a physical  separation  of  less  than  one  hundred 
kilometers  at  Transit  altitude  (1100  kilometers).  Two  satellites  are  in 
the  same  circular  orbit  while  the  third  orbit  is  slightly  eccentric. 

A.  Parameter  Reduction 

(i)  Covariance  Analysis 

Orbital  solutions  based  on  both  one  and  two  days  of 
synthetic  range  difference  data  were  generated  to  determine  what  effects 
multiple  satellite  data  processing  would  have  on  the  paramet.fr 
covariance.  The  data  noise  level  was  twenty  centimeters.  For 
individual  ortLt  processing  using  one  day's  data  from  five  stations  the 
one  sigma  orbit  parameter  uncertainty  implies  an  orbit  error  of  less 
than  one  meter.  The  one  sigma  uncertainty  in  the  determination  of  a 
drag  coefficient  over  this  span  of  data,  however,  implies  an  orbit 
error  of  about  three  meters.  Simultaneous  processing  in  which  a common 
drag  parameter  for  all  satellites  is  used  and  in  which  a single 
tropospheric  refraction  scaling  parameter  per  pass  is  allowed  yields  a 
one  sigma  uncertainty  in  orbit  parameters  of  some  ten  to  twenty  percent 
smaller  than  for  individual  processing.  The  one  sigma  drag  contribution 
to  orbit  error  was  reduced  to  1.7  meters.  These  results  indicate  that 
simultaneous  processing  is  capable  of  reducing  orbit  error  due  to  drag. 
This  same  analysis  was  done  for  two  days  of  data.  The  same  qualitative 
features  resulted.  It  should  be  emphasized  that  this  analysis  is 
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strictly  a covariance  analysis.  That  is  to  say  the  model  (geopotential, 
troposphere,  etc.)  is  assumed  perfect.  This  fact  is  the  origin  of  the 
very  small  uncertainties  quoted  above.  These  results  are  summarized  in 
Table  1. 

(ii)  Deterministic  Analysis 

The  effects  of  refraction  and  gravity  errors  on  both 
individual  and  simultaneous  orbit  determination  have  been  examined  over 
a two  day  period.  Orbit  solutions  were  based  on  synthetic  range 
difference  data  with  a standard  error  of  twenty  centimeters.  Orbit  so- 
lutions Indicate  that  tropospheric  refraction  bias  and  gravity  modeling 
(truncation)  error  which  is  recovered  primarily  through  the  drag 
parameter  generally  result  in  slightly  smaller  residuals  of  fit  for 
individual  processing  than  for  simultaneous  processing.  Individual 
satellite  processing  incorporating  both  a drag  parameter  per  satellite 
and  a refraction  bias  parameter  per  pass  is  better  able  to  "absorb" 
these  modeling  errors. 

However,  in  the  presence  of  these  unmodeled  errors 
simultaneous  processing  is  better  able  to  predict  the  relative  separation 
of  the  satellites  during  a three  day  period  following  the  two  day  fit 
span.  Table  2 gives  the  error  in  separation  distance  which  occurs  using 
both  processing  modes  for  one  pair  of  satellites.  As  indicated  the 
multiple  satellite  processing  scheme  predicts  the  relative  separation 
more  accurately.  This  is  graphically  depicted  for  the  case  of  trunca- 
tion of  the  spherical  harmonic  representation  of  the  gravity  field  to 
the  12,  12  level  in  Figure  1.  During  a three  day  prediction  the  relative 
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ORBIT  UNCERTAINTY  FOR  DIFFERENT  PROCESSING  METHODS 
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TABLE 


geometry  was  predicted  more  accurately  when  multiple  satellite  data 
processing  was  utilized. 

B.  Perturbation  Solution  Approximation 

This  method  was  first  applied  to  the  two  satellites  in  a 
common  circular  orbit.  A perturbation  solution  for  the  trailing 
satellite  (satellite  2)  was  developed  by  merging  its  integrated  inertial 
trajectory  with  the  integrated  perturbation  solution  from  the  leading 
satellite  (satellite  1).  This  was  accomplished  using  the  procedure 
described  above.  The  updated  position  for  satellite  1 at  time  t^ 
differed  by  less  than  ten  meters  from  the  updated  position  of  satellite 
2 at  time  t2«  Seven  orbit  fits  (one  iteration)  were  then  performed 
using  this  synthetic  perturbed  trajectory  to  determine  the  error  levels 
introduced  by  the  method.  These  solutions  were  based  on  range 
difference  observations  over  a two  day  period.  No  additional  modelling 
errors  were  introduced.  For  each  case  errors  in  the  data  were  intro- 
duced through  an  error  in  one  initial  orbit  constant  at  a level  compati- 
ble with  or  greater  than  the  poorest  initial  condition  errors  in  the 
precise  Transit  ephemeris  determined  at  NSWC  Dahlgren  Laboratory. 

Table  3 gives  the  maximum  orbit  errors  during  both  the  fit  span  and  a 
three  day  prediction  period  as  a function  of  the  orbit  parameter 
perturbed.  These  results  indicate  that  this  method  for  constructing  a 
synthetic  perturbed  trajectory  introduces  little  error  for  the 
satellites  in  a common  circular  orbit.  In  this  case  the  time  difference 
between  t^  and  t^  was  about  twelve  seconds. 
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TABLE 
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This  same  experiment  was  repeated  using  satellite  1 and 
satellite  3.  The  orbits  in  this  case  differ  to  a greater  extent  than 
the  first  pair  considered.  Using  the  same  technique  a synthetic 
perturbed  trajectory  was  created  for  satellite  3.  The  minimum  distance 
between  the  satellites  was  about  50  kilometers.  In  this  case  the  time 
differerce  between  t^  and  t^  was  about  60  seconds.  The  error  levels 
introduced  for  this  case  are  also  given  in  Table  3 . 
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Perturbation  Solution  Approximation  Error  Estimates 
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IV.  Sunmary  and  Recommendations 


This  report  has  examined  two  aspects  of  orbit  determination 
applicable  for  neighboring  satellites: 

(1)  The  simultaneous  processing  of  orbital  data  using  the 
reduced  parameter  set  approach  will  produce  orbits  with  an  accuracy 
consistent  with  individual  satellite  orbit  determination.  This  method 
however  yields  a better  solution  for  the  relative  separation  of  the 
satellites  in  the  presence  of  modeling  errors. 

(2)  The  use  of  perturbation  solutions  from  the  integration  of  one 
satellite's  variational  equations  to  process  other  satellites  has  been 
examined.  Although  this  approach  introduces  additional  error  into  orbit 
solutions  it  can  substantially  reduce  computation  time  required  to 
produce  these  orbits. 

Based  on  this  analysis  it  is  recommended  that  a multiple  satellite 
data  processing  scheme  be  adopted  for  satellites  in  neighboring  orbits. 
In  the  interest  of  system  accuracy  however  the  use  of  any  perturbation 
solution  approximations  should  be  limited  unless  computation  time 
requirements  become  critical. 
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Pass  Matrix  Concept 


Data  Class:  Range  Difference  (Integrated  Doppler) 

The  range  difference  data  class  is  given  by 


D( 


:.t  ) = P(t  ) - p(t  ) f 7 (tn-t  ) 
n n n-1  B ts  “tea 

+ 'f  7 (t_  - t ,) 2 + (1  + CL)  AR(t  ) - AR(t  )'  (A-la) 

B f n n 1 R n tea  j 


where 


P <tn>  = Ka^V’W 


(A- lb) 


and  where  the  vectors  rgat  and  rg  are  the  radius  vectors  of  the  satellite 
and  observing  station  respectively.  The  terms  fB  and  fg  represent  the 
frequency  bias  and  drift  characteristic  to  the  pass.  CR  represent  a 
refraction  scaling  parameter  and  AR  the  Hopfield  tropospheric  refraction 
correction  (Reference  2).  The  constant  fs  is  the  effective  frequency 
of  the  satellite,  c the  velocity  of  light,  and  tt  the  time  of 
closest  approach  of  the  satellite  to  the  tracking  station. 

Data  Aggregation 

The  pass  matrix  concept  (Reference  3)  as  used  in  orbit  determination 
involves  aggregating  satellite  tracking  data  on  a pass  basis.  Normal 
equations  are  formed  from  the  tracking  data  of  each  separate  pass  over  a 
station.  The  parameters  of  fit  are  the  satellite  orbital  constants  and 
certain  bias  parameters  characteristic  to  the  pass.  These  bias  para- 
meters represent  three  components  of  station  position,  refraction  bias, 
range  or  frequency  bias,  and  range  or  frequency  drift. 
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Formation  of  Pass  Normal  Equations  for  Satellite  "k"  and  Station  "A11 
After  data  from  the  j'th  pass  has  been  edited  to  eliminate  "bad" 
points,  the  normal  equations  are  formed: 

Let  Dt*  represent  the  data  taken  at  time  t^  with  associated  standard 
error  C£.  The  A matrix  Is  then  given  by 
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Table  A-l.  The  weight  matrix  % for  the  data  of  pass  j is  given  by 
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With  the  A^j  and  W 
for  pas  8 j and  given 


matrices  defined  the  least  squares  normal  equations 
by: 
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PARAMETERS  OF  FIT 


NOTATION 
Dynomicol  Terms : 

P, 

Pz 

P, 

P . 

P5 

Pe 

Cd 


PARAMETER 

x or  a 
y or  e sin  OJ 
i or  e cos  cj 
tot  I 

jf  or  L+G 

i or  Q 
Drag 


Bios  Terms : 
*.y,t 

Cr 

R g Or  fg 
Rg  Or  fg 


Station  Position 
Refraction  Bios 
Range  or  Frequency  Bios 

Range  or  Frequency  Drift 


TABLE  A- 1 
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where  the  nxl  vector X contains  the  observational  residuals.  The 

kj 


pair  [Bfcj,  T5kj]^  are  denoted  as  the  pass  matrix  for  the  j'th  pass. 
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Multiple  Satellite  Least  Squares 


Least  Squares  Method 


Let  [B^j,  denote  the  pass  matrix  (Appendix  A)  formed  from 

the  edited  observations  from  the  j'th  pass  of  satellite  k(k-l,2,3)  over 
station  l . These  normal  equations 


BkJ  iPW  ■ Ekj  •'.1,2.3 
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are  formed  for  the  orbital  and  bias  parameters  associated  with  the 
satellite-station  pair  k£  (Table  A-l).  A typical  pass  matrix  will  be 
of  the  form: 
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I Bb,b  ^b 


where  the  subscript  "o"  refers  to  the  orbital  parameters,  Cq  refers  to 
the  drag  scaling  parameters,  and  the  subscript  "b"  refers  to  pass  de- 
pendent bias  parameters. 
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For  each  pass  j the  elements  of  the  pass  matrix  (k-1,2,3) 

are  permuted  so  that  certain  pass  dependent  parameters  (station  and 
frequency  bias)  may  be  formally  eliminated.  The  equations  for  the 
elimination  of  these  pass  dependent  parameters  are 
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where  the  matrices  Ba,...,Bg  and  E^  and  are  constructed  by  sub- 
division of  the  matrices  of  equation  (B-2): 
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pair  kj  (k=l,2,3)  we  now  have  the  eliminated  pass  matrices  [B * * ]kjjt 

for  the  parameter  set  . 

Now,  under  the  assumption  that  the  atmospheric  drag  scaling  con- 
stants and  tropospheric  refraction  correction  parameters  per  pass  are 

common  among  all  three  satellites,  the  normal  equations 
for  the  j'th  pass  may  be 
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combined  to  yield  the  combined  pass  matrix  [B1,  E']^: 

CbM'], 


K,0 


ji 

2 


i 


i 0 ! 0 i b°,cd1 
r i • 

! B'2  ! 0 ! B'  2 

I o»°  ' ! °»CD 


l 


i 


I 


3 ‘ 

IB'  ! B'  3 

l o,o  ; o»Cjj 


7 ,k 

£-B'  £L 

I r Cd.^d 


’W  : 


E\ 


B' 


L. 


1 

2 1 
1 CR  1 

J 

E'2 

0 

3 1 

c*  i’ 

!■> 

O 

i ■ 

r- 1 

k i 

cd»cr  1 

i 

L T 1 

k CD 

;•  k 1 

Si’S  I 

k °R 

At  this  point  the  elimination  equations  may  again  b'e  utilized  to 
eliminate  the  pass  dependent  tropospheric  refraction  scaling  parameter 
Cr  yielding  a new  set  of  eliminated  pass  matrices  involving  only 
dynamic  parameters: 

BV  EP"  - E"  . (B-7) 

J (19x19)  (19x1)  3(19x1) 


At  this  point  all  pass  normal  equations  (B-7)  may  be  summed  to  yield  the 
final  arc  normal  equations  for  the  solution  of  dynamical  parameters 
(six  per  satellite  along  with  a common  drag  for  the  three  satellites) : 
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